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ABSTRACT 

Context. While planets in the solar system only have a low inclination with respect to the ecliptic there is mounting evidence that 
in extrasolar systems the inclination can be very high, at least for close-in planets. One process to alter the inclination of a planet is 
through planet-disc interactions. Recent simulations considering radiative transport have shown that the evolution of migration and 
eccentricity can strongly depend on the thermodynamic state of the disc. So far, this process has only been studied for a few selected 
planet masses using isothermal discs. 

Aims. We extend previous studies to investigate the planet-disc interactions of fixed and moving planets on inclined and eccentric 
orbits. We also analyse the efl'ect of the disc's thermodynamic properties on the orbital evolution of embedded planets in detail. 
Methods. The protoplanetary disc is modelled as a viscous gas where the internally produced dissipation is transported by radiation. 
To solve the equations we use an explicit three-dimensional (3D) hydrodynamical code NIRVANA that includes full tensor viscosity, 
as well as implicit radiation transport in the flux-limited diffusion approximation. To speed up the simulations we apply the FARGO- 
algorithm in a 3D context. 

Results. For locally isothermal discs, we confirm previous results and find inclination damping and inward migration for planetary 
cores. For low inclinations (ii:2H/r), the damping is exponential, while it follows di/dt oc t- for larger /. For radiative discs, the 
planetary migration is very limited, as long as their inclination exceeds a certain threshold. If the inclination is damped below this 
threshold, planetary cores with a mass up to a 33Me^i„i, start to migrate outwards, while larger cores migrate inwards right from the 
start. The inclination is damped for all analysed planet masses. 

Conclusions. In a viscous disc an initial inclination of embedded planets will be damped for all planet masses. This damping occurs 
on timescales that are shorter than the migration time. If the inclination lies beneath a certain threshold, the outward migration in 
radiative discs is not handicapped. However, only planets with a mass up to « 33Me„„i, are prone to this outward migration. Outward 
migration is strongest for circular and non-inclined orbits. 

Key words, accretion discs - planet formation - hydrodynamics - radiative transport - planet disc interactions - inclination 



1. Introduction 

Planets in our solar system do not all move in the same plane, 
but are inclined with respect to the ecliptic. The inclinations 
are in general small. The highest inclined planet is Mercury 
(i - 7.01°), and only some of the dwarf planets are inclined 
much higher (e.g. Pluto at ; = 17°, or Eris at / - 44.2°). This 
low inclination with respect to the ecliptic is typically taken as 
an indication that planets form within a flattened protoplanetary 
disc. The observed high / for smaller objects are then a result of 
gravitational scattering processes in the early solar system evo- 
lution. Recently, it has been discovered that a relatively large 
number of transiting close-in planets in extrasolar systems can 
be inclin ed significantly wi th respect to the stellar rotation axis 
(see e.g. iTriaud et al.ll2010l) . This discovery makes it necessary 
to take inclination into consideration in both numerical simula- 
tions and theoretical consideratio ns of planet-disc interac t ion. 

Starting from early work (Gol dreich & Tremaine ] 119801 

iPapa loizou & Lin 1984; Ward 1986), the effect a protoplanetary 
disc has on the migration of the planets has been analysed in- 
tensively aft er the discovery of hot J upiters close to the star (for 
a review see IPapaloizou et aTl (l2007l) ). Numerical simulations of 
planet-disc interactions are a valuable tool for a better under- 
standing. Two-dimensional ( 2D) s tudies have been perfor med 
for example bv lBrvden et"al1(ll999b : iD'Angelo et"ai1 (|2002|) and 
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later f ully three d imensional (3D) simulations (ID'Angelo et al.l 
200l iBate et al.l 12003) followed. All these simulations only 
considered isothermal discs. However, recent s ttidies (start- 
ing with the work of iPaardekooper & Mellemal (l2006l) ) have 
shown that including radiation transport in planet-disc interac- 
tion studies resulted in a slowed down or even reversed mi- 
gratio n iBaaiteau & Masset 2008; Paardekooper & Papaloi zou! 
2008; Paardekoo per & Mellemal 120081; iKlev & Cridal 120081 
Avliff'e & BatQ^OloT All authors agree that the inclusion of ra- 
diative transfer is important and that it strongly affects the migra- 
tion properties of low mass planets (type-I-migration) in discs, 
even though there are differences in the magnitude of the effect. 
It also affects whether it actually leads to reversal of migration 
or only to a slow down. 

Recently, we have analysed the dis c's influence on the or- 
bital eccentricity (iBitsch & Kleyl l20Tol) for planets of different 
masses in fully radiative discs. We showed that planet-disc inter- 
action leads quite generally (for all planet masses) to a damping 
of planetary eccentricity on shorter timescales than the migra- 
tion timescale. Additionally, planets on eccentric orbits embed- 
ded in fully radiative discs only migrate outwards, when the ec- 
centricity is damped to a nearly circular orbit. Clearly, eccentric- 
ity destroys the very sensitive corotation torques near the planet, 
which results in inward migration until the orbit is circularised 
again. 

The effect of disc-planet interaction on inclination has not 
yet been studied in any wide extent, owing to the heavy compu- 



1 



Bitsch & Kley: Evolution of inclined planets in three-dimensional radiative discs 



tational requirements. iTanaka & Wardl(l2004l) analysed the influ- 
ence of the disc on the inclination of low mass planets in linear 
studies. They find an exponential damping for any non-vanishing 
inclination, but their results are only formally valid for ; « H/r. 
Numerical simulations of inclined planets with higher inclina- 
tions than H/ r (inclination in radians) have shown that exponen- 
tial damping may be valid up to / ^ 2H/ r. For planets on higher 
initial inclined orbits, the damping rate deviates from being ex- 
ponential, and it can b e fitted best with a dildt oc function 
(iCresswell et al ] |2007h . Only low mass planets with a mass of 
about 20-30 earth masses were considered in these studies. The 
in fluence on an inc lined Jupiter type planet has been considered 
bylM arzari & Nelson (2009.). They find that highly inclined and 
eccentric planets with Jovian masses lose their inclination and 
eccentricity very fast (on a timescale of order of 10^ years), when 
entering the disc again. Since a highly inclined planet is only dis- 
turbed in a small way by the accretion disc (and vice versa), such 
a planet is only able to op en a gap in the disc , when the inclina- 
tion drops to ; < 10.0°. iTerquem & Aimial (1201 Ol) considered 
inclination evolution utilising a Kozai-type of effect for planets 
on high inclined and eccentric orbits in the presence of the disc. 

In the present work we will extend the work on inclined plan- 
ets and analyse the effect of planet-disc interaction on embedded 
planets of all masses. In fully three-dimensional hydrodynamical 
simulations we will consider isothermal as well as fully radiative 
discs and study the change in inclination of an embedded planet. 

In Section |2] we give a short overview of our code and nu- 
merical methods. In Section|3]we measure the change of inclina- 
tion and semi-major axis for IQMEanh planets on fixed circular 
inclined orbits and then let these planets evolve their orbits in 
the disc due to the torque acting on the planet form the disc in 
Section H) In Section |5] we study the influence of the planetary 
mass on the evolution of planets on inclined orbits. In Section 
|6] we follow the evolution of IQMEanh planets on eccentric and 
inclined orbits. In Section |3] to Section |6] we also point out the 
differences between the isothermal and fully radiative simula- 
tions. Finally we summarise and conclude in Section|7] 



2. Physical modelling 

The protoplanetary disc is modelled as a three-dimensional 
(3D), non-self-gravitating gas whose motion is described by 
the Navier-Stokes equations. We treat the disc as a vis- 
cous medium, where the dissipative eff'ects can then be de- 
scribed via the standard viscou s stress-tensor approach (eg. 
iMihalas & Weibel MihalaslfTgsl . We also assume that the heat- 
ing of the disc occurs solely through internal viscous dissipa- 
tion and ignore the influence of additional energy sources (e.g. 
irradiation form the central star). This internally produced en- 
ergy is then radiatively diffused through the disc and eventually 
emitted from its surface. For this process we use the flux-limited 
diffusion approximation (FLD, Levermore & Pomraning 1981), 
which allows us to treat the transition from optically thick to 
thin regions as an approximation. A more detailed description of 
the modelling an d the numer i cal me thodo logy is provided in ou r 
previous papers (iKlev et al.l [20091) and (iBitsch & Kiev! l201(]h . 
Compared to our previous papers, we now extend the simula- 
tions by including planets with different masses on inclined or- 
bits. 



2. 1 . General Setup 

An important issue in modelling planetary dynamics in discs is 
the gravitational potential of the planet since this has to be arti- 
ficially smoothed to avoid singularities. 

While in two dimensions a potential smoothing takes care of 
the otherwise neglected vertical extension of the disc, in three 
dimensional simulations the most accurate potential should be 
used. As the planetary radius is much smaller than a typical grid 
cell, and the planet is treated as a point mass, a smoot hing of the 
potent ial is required to ensure numerical stability. In iKlev et all 
(l2009h we have discussed two alternative prescriptions for the 
planetary potential. The first is the classic e-potential 



Of, = - 



Gm„ 



(1) 



Here mp is the planetary mass, and d -\r - rp| denotes the dis- 
tance of the disc element to the planet. This potential has the 
advantage that it leads to very stable evolutions when the pa- 
rameter e (stated in units of the Hill radius) is not too small, i.e. 
is a significant fraction of the Roche-lobe size. The disadvan- 
tage is that for smaller e, which would yield a higher accuracy 
at larger d, the potential becomes very deep at the planetary po- 
sition. Additionally, due to the long range nature, the potential 
differs from the exact one even for medium to larger distances d 
from the planet. 

To resolve these problems at small and large d simul- 
taneously, we have suggested the following cM^j/c-potential 
(iKlahr & Klevll200^ : lKlev et alJl2009l) 
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Here, rj„, is the smoothing length of the potential in units of the 
Hill radius. The construction of the planetary potential is in such 
a way that for distances larger than r „„ the potential matches the 
correct 1/r potential. Inside this radius {d < rsm) it is smoothed 
by a cubic polynomial. This potential has the advantage of exact- 
ness outside the specified distance r„„, while being finite inside. 
The parameter r„„ is equal to 0.5 in all our simulations, unless 
stated otherwise. The value of the smoothing length of the plan- 
etary p otential Fs,,, was discussed in great detail in iKlev et all 
(|2009|) . where we compared different smoothing lengths and po- 
tential types, resulting in our choice of r^,,, = 0.5. In this work we 
use the cubic-form for the planetary potential for all fully radia- 
tive simulations and for the isothermal simulations with planets 
on fixed orbits. We note that for e - 0.25 the depth of the e- 
potential at the position of the planet, d = 0, is identical to the 
cubic -potential with r^,,, = 0.5. However, the difference reaches 
a maximum of about 20% at li = 0.25, is 10% at li = 0.5, and 
slowly diminishes for larger d. 

Our simulations for moving planets in isothermal discs have 
shown that in this case the r,,,, = 0.5 cubic potential yields 
unstable evolutions. This is due to the relatively steep centre 
which leads to a significant density enhancement for isothermal 
discs. Hence, we used in this case directly the much smoother 
e-potential with e = 0.8 instead, which has been used before in a 
similar context dCresswell et al.ll2b07h . without having explored 
other choices here. 

The gravitational torques and forces acting on the planet are 
calculated by integrating over the whole disc, where we apply a 
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tapering function to exclude the inner parts of the Hill sphere of 
the planet. Specifically, we use the smooth (Fermi-type) function 



Md) = 



dIRH-b 



(3) 



which increases from at the planet location {d = 0) to 1 out- 
side d > Rh with a midpoint ft - 1/2 at (i = bRH, i.e. the 
quantity b denotes the torque-cutoff radius in units of the Hill 
radius. This torque-cutoff is necessary to avoid large, probably 
noisy contributions from the inner parts of the Roche lobe and 
to disregard material that is possibly gravitationally bound to the 
planet (jCrida et al. 2009i). Here we assume (as in our previous 
paper) a transition radius of b - 0.8 Hill radii, as a chang e in b 
did not influence the results significantly ( iKlev et alj|2009t) . 

2.2. Initial Setup 

The three-dimensional (r, 6, <p) computational domain consists 
of a complete annulus of the protoplanetary disc centred on 
the star, extending from r„,„, = 0.4 to r„,ax - 2.5 in units of 
'"o = <ijup = 5.2AU. In the vertical direction the annulus extends 
7° below and above the disc's midplane, meaning 83° < 9 < 97°. 
Here 6 denotes the polar angle of our spherical polar coordinate 
system measured from the polar axis. The central star has one 
solar mass M* - Mq, and the total disc mass inside [r„„„, r,„av] 
is Mdisc = O.OIMq. For the isothermal simulations we assume 
an aspect ratio of H/r - 0.037 for the disc, in very close agree- 
ment with the fully radiative models of our previous studies. For 
the radiative models H/ r is calculated self-consistently from the 
equilibrium structure given by the viscous internal heating and 
radiative diff'usion. The isothermal models are initiahsed with 
constant temperatures on cylinders with a profile T{s) oc s^' 
with s = r sin 9. This yields a constant ratio of the disc's vertical 
height H to the radius s. The initial vertical density stratification 
is approximately given by a Gaussian: 



p(r, 9) = po(r) exp 



(7r/2 - of 



(4) 



Here, the density in the midplane is po{r) oc r which leads 
to a 2(r) oc r"'^^ profile of the vertically integrated surface den- 
sity. In the radial and ^-direction we set the initial velocities to 
zero, while for the azimuthal component the initial velocity is 
given by the equilibrium of gravity, centrifugal acceleration and 
the radial pressure gradient. This corresponds to the equilibrium 
configuration for a purely isothermal disc. 

For our fully radiative model we first run a 2D axisymmet- 
ric model (starting from the given isothermal equilibrium) to ob- 
tain a new self-consistent equilibrium where viscous heating bal- 
ances radiative transport/cooling from the surfaces. After reach- 
ing that equilibrium, we extend this model to a full 3D simula- 
tion, by expanding the grid into ^-direction. The resulting disc 
for this model has H/r ^ 0.037 so we choose that value for our 
isothermal runs. 



2.3. Numerical Setup 

Our coordinate system rotates at the initial orbital frequency of 
the planet (at r = ro). We use an equidistant grid in r, 9, (p with 
a resolution of (A^^, Ne, N^) = (266, 64, 768) active cells for our 
simulations. At r„,/„ and r,„ax we use damping boundary condi- 
tions for all three velocity components to minimise disturbances 
(wave reflections) at these boundaries. The velocities are relaxed 



towards their initial state on a timescale of approximately the lo- 
cal orbital period. The angular velocity is relaxed towards the 
Keplerian values, while the radial velocities at the inner and 
outer boundaries vanish. Reflecting boundary conditions are ap- 
plied for the density and temperature in the radial directions. We 
apply periodic boundary conditions for all variables in the az- 
imuthal direction. In the vertical direction we set outflow bound- 
ary conditions for 0,„„, = 83° and 0„,at = 97° (t he surfaces of the 
disc). We use the finite volume code NIRVANA (IZiegler & Yorkd 
Il997h with implicit radiative transp ort in the flux- hmited diffu- 
sion approx imation and the F ARGO (lMassetll2000l) extension as 
described in lKlev et al.l(l2009l) . 

2.4. Simulation Setup 

In the first part of our model sequence we consider the or- 
bital evolution of a planet with a fixed mass (IQMEarth) on in- 
clined orbits using different initial inclinations. For comparison 
we consider isothermal and fully radiative models. Using radia- 
tive discs here is a direct extension of a previous study under 
purely isothermal disc conditions using the same planet mass 
dCresswell et al.ll2007h . We distinguish two different approaches 
for these lOMEanh models: first, a model sequence where the 
planet stays on a fixed inclined orbit and secondly where the 
planet is free to move inside the computational domain under 
the action of the planet-disc gravitational forces. For the sec- 
ond models we insert the planet in the disc and let it move 
immediately, but using a time-dependent mass growth of the 
planet (through the planetary potential) with a typical switch 
on time of 10 planetary orbits. For the first set of models the 
IQMEanh planet is inserted at the nominal mass in the disc at 
the start of the simulation. Initially the planet starts at a dis- 
tance r = ajup = 5.2AU from the central star. For the fully 
radiative simulations we set the ambient temperature to a fixed 
value of IQK at the disc surface (at 9„un and 9,„a.^), which en- 
sures that all the internally generated energy is liberated freely 
at the disc's surface. This low temperature boundary condition 
works very well at optically thin boundarie s and does not ef- 
fect the inner par ts of the optically thick disc (iKlev & Linll 19991 : 
iKlev et aT]|2009h . In the second part of the project we consider 
sequences of models for a variety of planet masses. We note that 
a lOMEarth planct has in our simulations, using our standard res- 
olution, a Roche radius of about 3.3 grid cells. In the last part 
we limit ourselves to a lOMEanh planet and simulate the evolu- 
tion of planets with an initial eccentricity and inclination. For all 
isothermal simulations, which allow the planet to move freely 
inside the disc, we use the e-potential. Only for isothermal sim- 
ulations of planets on fixed inclined orbits, the cubic potential is 
used. Additionally, the cubic potential is used for all fully radia- 
tive simulations as well. 



3. Models with an embedded planet on fixed 
circular inclined orbits 

In this section we consider planets on fixed circular and in- 
clined orbits embedded in either isothermal or fully radiative 
discs. From the disc forces acting on the planet on the fixed or- 
bit we calculate its change of inclination and its migration rate. 
In this section all simulations use the cubic potential, featuring 
Mpianei = ^.QMEariiu and a scmi-major axis of a = 1.0. Our simu- 
lations only cover up to 7° above and below the equatorial plane, 
as the disc gets very thin for these regions. However, we investi- 
gate the motion of planets for higher inclinations (for fixed plan- 



3 



Bitsch & Kley: Evolution of inclined planets in three-dimensional radiative discs 



ets up to 15°). The results are not influenced by our limited ver- 
tical extent of our computational grid, as the density is very low 
in the upper layers of the disc. 

3.1. Change of inclination 

To determine the change o f orbita l elements for planets on fixed 
inchned orbits, we follow iBurns I (119761) . If a small disturbing 
force dF (given per unit mass) due to the disc is acting on the 
planet, the planet changes its orbit. This small disturbing force 
dF may change the planetary orbit in size (semi-major axis a), 
eccentricity e and inclination /. The inclination ; gives the angle 
between the orbital plane and some arbitrary fixed plane, which 
is in our case the equatorial {6 = 90°) plane, which corresponds 
to the midplane of the disc. Only forces lying in the orbit plane 
can change the orbits size and shape, while th ese forces can not 
change the orientation of the orbital plane. In iBurnsI (1 19761) the 
specific disturbing force is written as 

dF = R + T -H N = iJcR + rer + A^e^v , (5) 

where the e's represent an orthogonal unit vector triad. The per- 
turbing force can be split in these components: R is radially out- 
wards along r, T is transverse to the radial vector in the orbit 
plane (positive in the direction of motion of the planet), and N is 
normal to the orbit planet in the direction R x T. 
[Burn s ( 1976) finds for the change of incUnation 

di aN cos ^ 



H 



(6) 



where the numerator is the component of the torque which ro- 
tates the specific angular momentum H = r x f about the line of 
nodes (and which thereby moves the orbit plane). The specific 
angular momentum H is defined as 

H = yjGMail . (7) 

The angle ^ is related to the true anomaly / in the following 
way f - ^ - u, with to being the argument of periapsis. ^ de- 
scribes the angle between the line of nodes and the planet on its 
orbit around the star For our case of circular orbits the argument 
of periapsis u is zero. 

In the simulations the planet's inclination spans from / - 
0.5° to / = 15.0°, in the isothermal as well as in the fully radia- 
tive regime. In Fig.[T]we display the rate of change of inclination 
for planets on fixed inclined orbits in the isothermal and fully ra- 
diative scheme, after we average dildt over 2 planetary orbits, 
after 150 orbits. As the angle ^ changes after every time step, 
one needs to average the quantity of dildt over the time of a 
planetary orbit to determine an exact value for the change of in- 
clination of the planetary orbit. For both thermodynamic systems 
the change in inclination is nearly identical and always negative, 
meaning that moving planets with a non zero inclination lose 
their initial inclination at a rate according to the current inclina- 
tion. At / = 4.0° the loss of inclination is at it's maximum and 
is reduced for higher and lower inclinations. For high inclined 
planets (/ > 6.0°) this loss is difdt oc as indicated by the 
short-dashed fit in black in Fig.[Tl and for lower inclined planets 
(i < 4.0°) the inclination is damped with an exponential decay as 
indicated by the linear fits in light bl ue and purple in F'ig. [T] This 
is identical to the behaviour found in lCresswell et al.l(l2007h and 
is also confirmed by our simulations of moving planets below. 

From our measured di/dt for planets on fixed orbits, we can 
calculate the inclination damping timescale t,„c according to 
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Fig. 1. Calculated rate of change of the inclination (^) for in- 
clined lOMEarth plancts on fixed circular orbits. Results for the 
isothermal (red) and fully radiative simulations (dark blue) are 
shown. Overlayed are in light blue and purple fits for the expo- 
nential decay of inclination for small inclinations and in black 
for the decay di/dt oc for larger inclinations. 



with ; being the inclination of the planet. In the exponential, 
small / regime we obtain for the isothermal simulations t,„(. = 27 
orbits and for the fully radiative simulations T „e = 20 orbits. 
Below we compare this to the linear results of iTanaka & WardI 
(l200l and moving planets. 



3.2. Ctiange of semi-major axis 

The power on the planet determines the change in semi-major 
axis of the planet, while the torqi ie represents a chang e in both 
eccentricity and semi-major axis (iBitsch & Klevlf2010l) . For cir- 
cular orbits torque and power are identical. For that reason we 
display and refer only to the torque acting on the planet, instead 
of mentioning the power as well. A positive torque indicates out- 
ward migration, while a negative torque represents inward mi- 
gration. 

In Fig.|2]we display the torque acting on the planet for both 
thermodynamic systems. The torque is positive for the fully ra- 
diative disc and negative for the isothermal disc. In the fully ra- 
diative scheme the torque has its maximum at ; = 0.0° with a 
massive loss towards higher inclinations until about ; » 4.5°, 
when the torque is about zero. For even higher inclinations the 
torque remains at about zero level. In the isothermal case the 
torque is negative and nearly constant for all inclinations, but 
with an increase towards higher inclinations. 

As for circular orbits the torque determines the change of 
semi-major axis for our embedded planets, we expect a positive 
migration rate for planets in a fully radiative disc and a negative 
migration rate for planets in an isothermal disc. To det ermine 
a mig ration rate for planets on fixed orbits, we follow iBumsl 
(1 19761) again: 



da 



dt {GM,yi^ 



.3/2 



(\-e^y^l^[Resmf + T{\+esinf)\ , (9) 



dildt 



1 



(8) 



where a is the semi-major axis, e the eccentricity, f - ^ - to the 
true anomaly and R and T are radial and tangential directions 
of the disturbing force, see eq. (|5). Please note that only forces 
lying in the orbit plane can change the orbit size, and that for 
circular orbits the specific total torque on the planet is simply 
given by Ytot - aT, with H = r,o,. 
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Fig. 2. Torque acting on a lOMEanh planet on a fixed inclined cir- 
cular orbit in dependency of the inclination and thermodynamics 
of the disc. Torque and power are the same for planets on circu- 
lar orbits. The torque has been averaged over 20 orbits, taken 
from t - 180 to f = 200 orbits for the isothermal (solid red line) 
and from t = 120 to f - 140 for the fully radiative simulation 
(dashed blue hne). 



The rate of migration is displayed in Fig. |3] In the fully ra- 
diative scheme the planet experiences a positive migration rate 
(outward migration) for all inclinations smaller than / < 4.5°, but 
with a strong increase towards lower inclinations (about a factor 
of 10 for the difference between ; = 4.0° and / = 1.5°). For 
inclinations lower than / x 1.5° the migration rate stays approx- 
imately constant. Planets with an inclination of 4.5° and higher 
seem to be stalled and migrate neither inwards nor outwards. The 
effects responsible for outward migration seem to be only valid 
in the small inclined case. 

In the isothermal case the migration rate is nearly constant 
and negative with a slight increase towards higher inclinations, 
suggesting that an inclined planet will move inwards in the 
isothermal regime for low inclined planets. Higher inclined plan- 
ets will move inwards only at a very small rate compared to their 
low inclined counterparts. For large / the change of semi-major 
axis approaches zero, as the planet's orbit is high above the dense 
regions of the disc, which are capable of influencing the planet's 
orbit. 

From these results we conclude that a planet with a non-zero 
initial inclination will lose this inclination in time. This inclina- 
tion loss has in principal no effect on the trend of migration, so 
a planet in a fully radiative disc will migrate outwards, while a 
planet in an isothermal disc will migrate inwards. In one of the 
next chapters we will observe moving planets with a non-zero 
inclination, which will do exactly that. 



3.3. Torque analysis 

To understand the behaviour of the total torque in more detail 
we now analyse the space-time variation of the torque of the 
planet. For that purpose we introduce the radial torque density 
r(r), which is defined in such a way that the total torque r,oi 
acting on the planet is given by 



Fir) dr. 



(10) 
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The radial torque density has been a very useful tool to inves- 
tigate the origin of the torques in our previous work on planets 



2 4 6 8 10 12 14 
Inclination 

Fig. 3. Calculated rate of change of the semi-major axis (d/a) 
for inclined lOMEanh planets on fixed, circular orbits. For incli- 
nations lower than / < 4.5° we observe a positive migration rate, 
which suggest that low inclinations do not stop outward migra- 
tion, while we observe a negative migration rate for the isother- 
mal case. 



on fixed circular orbits. Even though the planet may be on an in- 
clined orbit, this does not alter our definition of the radial torque 
density r(r) which is calculated in any case with respect to the 
orbital plane of the planet. In Fig.|4|we display r(r) for a selec- 
tion of our planets in the isothermal and fully radiative regime. 

In the isothermal case the r(r)-function follows the same 
trend for all displayed inclinations. The only difference is near 
the planet, around the r = 1.0 region. The torque of the low 
inclined planets declines constantly, while the higher inclined 
planets (2.0 < / < 4.0) form a little spike in the r(r) distribu- 
tion, which disappears for even higher inclinations again. For 
higher inclinations the Lindblad-Torques become less and less 
pronounced, as can be seen in Fig.|4|(top panel). Please keep in 
mind that for the inclined planets the torque density changes dur- 
ing one orbit, as can be seen from the change in the total torque 
acting on the planet (top panel in Fig.|6|. This is also illustrated 
in Fig.|5]for the planet with / = 2.0°. It can been seen that for 
symmetry reasons r(r) is identical at upper and lower culmina- 
tion, and at the crossings of the nodal line, respectively. 

In Fig. |6|the torque acting on planets (please note that for 
circular orbits torque and power are identical) with different in- 
clinations on fixed circular orbits (top) and the normal compo- 
nent of the disturbing force dF (bottom) for fully radiative discs 
is displayed. The torque acting on the planets oscillates with 
time, making two oscillations every orbit. The planet starts at 
the highest point in orbit (upper culmination) and evolves to the 
lowest point in orbit (lower culmination) in the time of half an 
orbit. Then the planet moves to the upper culmination again and 
reaches it at the end of the orbit. The planets with a lower inclina- 
tion (/ < 4.0°) have clearly a positive average torque, while for 
the planets with higher inclinations positive and negative con- 
tributions approximately cancel out such that the average total 
torque is very small (see Fig. |2). Higher inclinations also trigger 
higher amplitudes in the torque distributions. These amplitudes 
in the torque distribution suggest oscillations in the evolution of 
the semi-major axis of the planet (to a small degree), which are 
visible when the planet is allowed to move freely inside the disc 
(see Fig.ini Fig.[T3]and Fig.[T5]l. 

In the lower panel in Fig. |6]the normal component of the dis- 
turbing force cos ^ is displayed. cos ^ oscillates also twice 
in every orbit, as the torque, but it is slightly shifted with respect 
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Fig. 4. Torque density acting on a planet on a fixed inclined or- 
bit in dependency of the inclination of the planet embedded in an 
isothermal Hjr - 0.037 disc (top) and a fully radiative disc (bot- 
tom). The snapshots are taken at f = 90 orbits for the isothermal 
simulations and at f = 150 orbits for the fully radiative simu- 
lations. At this time planet- disc interactions are in equilibrium 
and the planet is at its lowest point in orbit (lower culmination). 



to the torque. Negative values indicate a reduction of inclina- 
tion, while a positive force indicates an increase of inclination. 
The oscillations in cos ^ indicate oscillations in the inclination, 
when the planet is allowed to move freely in the disc. Indeed, in 
case the planet moves inside the disc, these oscillations in the 
inclination become visible at the beginning of the simulations 
(see Fig.[TTl Fig. [T3] and Fig.[T4lbelow). On average the normal 
component of the disturbing force cos ^ is negative, which in- 
dicates inclination damping for all simulated inclinations, see 

Fig.m 

For simulations in the isothermal regime, the results are 
nearly identical, with only one exception: The overall torque is 
negative, which results in inward migration. The normal compo- 
nent of the disturbing force, cos ^, is also very similar com- 
pared to the results in the fully radiative regime. We therefore 
forgo to discuss the isothermal results in detail at this point. 

To illustrate the dynamical impact of the planet on the disc, 
we display the surface density of planets with / = 1 .0° and ; - 
4.0° in Fig. in 

The surface density of the / — 1 .0° planet is very similar to 
that of a non inclined planet (not displayed he re, but compared 
to bottom picture of Fig. 8 in lKlev et alJ (l2009l) ). and the torque 
density r(r) is nearly identical. The surface density also shows 
no disturbances for the / = 1 .0° planet, so that we only observe 
the Lindblad Torques in the isothermal case, see Fig. |4] However, 




Fig. 5. Torque density acting on a 20M(.(„,/, planet on a fixed in- 
clined orbit with / = 2.0° in a fully radiative disc at different 
phases within its orbit. Integer (half integer) values of the indi- 
cated time refer to the upper (lower) culmination. At the inter- 
mediate quarterly orbit times the planet moves through the nodal 
line. 



the surface density of the higher inclined planet (/ = 4.0°) shows 
density disturbances near the planet. For the ; - 1 .0° we also ob- 
serve a higher density in the planets surroundings and a stronger 
pronounced spiral wave compared to higher inclined planets. 
The less pronounced spiral wave for the higher inclined planets 
leads to a smaller Lindblad Torque for these planets compared 
to the low inclined planets, which can be seen clearly in Fig.|4] 

In Fig. [8] the density of a r-0 slice through the planets in the 
isothermal Hjr - 0.037 disc is display. The planets are at their 
lowest point in orbit when the snapshot was taken. On the one 
hand, the ; = 1 .0° planet (top) has accumulated more mass in its 
vicinity as the ; = 4.0° planet (bottom), as the density in the mid- 
dle of the disc is higher and thus it is easier to accumulate mass. 
On the other hand, the higher inclined planet seems to disturb 
the density structure of the disc stronger than the lower inclined 
planet. This distortion in the density distribution in the r-0 plane 
reflects in the distorted surface density structure, displayed in the 
bottom picture of Fig. |7] 

In the fully radiative case, the torque density (Fig.|4) shows 
the well discussed spike in the torque distribution at r = 0.984 
(see Kiev et al. (2009)). This spike is more pronounced for lower 
inclinations (; < 2.0°) and is reduced for the higher inclinations. 
This reduction for the torque density at r = 0.984 causes the total 
torque to decrease for higher inclinations, see Fig. |2] For higher 
inclinations the Lindblad Torques are also reduced compared to 
the lower inclined planets. 

In Fig. |9] the surface density of IQMEarth planets in a fully 
radiative disc with / - 1 .0° and / = 4.0° are displayed. The over- 
all surface density distribution of the ; = 1.0° planet is nearly 
identical with the su rface densi t y dist ribution of a non inclined 
planet (see Fig. 12 in lKlev et alj (l2009h ). which is also supported 
by the identical torque density distributions. The surface density 
for the higher inclined planet (/ - 4.0°) on the other hand shows 
some differences. The density ahead of the planet (r < 1.0 and 
<p > 180°) is reduced and the density behind the planet (r > 1.0 
and (f) < 180°) is increased compared to the density distribution 
of the low inclined (/ - 1.0°) planet. This loss and gain of den- 
sity results directly in a reduced spike (at r - 0.984) in the torque 
density (Fig. |4]i, which reduces the total torque. The density in- 
side the planet's Roche lobe is also reduced in the high inclined 
case, as the planet can not accumulate as much mass as a lower 
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Fig. 6. In the top panel we display the evolution of the torque act- 
ing on planets in a fully radiative disc on fixed circular, inclined 
orbits in time. In the bottom panel the normal component of the 
disturbing force dF, cos ^ as to be used eq. is displayed in 
dependency of time. 



inclined planet. This feature was also visible in the isothermal 
case, but keep in mind that the planet in the isothermal case ac- 
cumulates more mass than in the fully radiative disc. Also the 
density in the spiral waves is reduced in the high inclined case, 
which results in lower Lindblad Torques compared to the low 
inclined case, see Fig.H] 

In Fig.[TO]the density through the fully radiative disc (a r-9 
slice) at the planets location for / = 1.0° and ; = 4.0° is dis- 
played. Clearly the higher inclined planet accumulated less mass 
than the lower inclined planet, as the density in the surround- 
ings of the high inclined planet is reduced compared to the discs 
midplane. Also the higher inclined planet seems to disturb the 
disc's density more than the lower inclined planet, which results 
in small fluctuations in the surface density, see Fig. |9] The lit- 
tle fluctuations in the surface density are less compared to the 
isothermal case, as the radiative transport/cooling smoothes the 
density in the disc and the planet does not accumulate as much 
mass as in the isothermal case, which disturbs the density struc- 
ture as well. 

In the r - slices of the disc, the difference in the mass ac- 
cumulation of the planet in the isothermal and fully radiative 
regime is very obvious. This effect is caused by the thermody- 
namics of the disc, as the isothermal disc is unable to heat and 
cool, the planet can accumulate more mass as in the fully radia- 




Fig. 7. Surface density for lOMEanh planets on fixed circular in- 
clined orbits embedded in an isothermal Hjr - 0.037 disc with 
/ = 1.0° (top) and / = 4.0° (bottom). The planet is at its lowest 
point in orbit (lower culmination) at the time of the snapshot. 



tive regime. Nevertheless, only in the fully radiative regime the 
calculated change of semi-major axis predicts outward migration 
for planets on fixed orbits with /q < 4.5°. 



4. Moving planets on initially inclined orbits 

In the previous section we calculated the change of inclination 
and semi-major axis for IQMEanh planets on fixed inclined orbits 
in the isothermal and fully radiative regime. The results stated 
that the planets will lose their inclination in time and will mi- 
grate inwards in the isothermal and outwards in the fully ra- 
diative scheme (for /q < 4.5°). We now want to confirm these 
results by allowing the planets to move freely inside the disc. 
For the following simulations we use our standard resolution, 
with a 2QM Earth planet embedded at r = 1.0 with different in- 
clinations ;. The discs thickness is H/r = 0.037. In the first 10 
orbits of the planetary evolution the mass of the planet will rise 
until it reaches its desired mass of 2QM Earth after 10 planetary 
orbits. This way the disc will not be disturbed as much as by 
putting a planet with its full mass inside the disc at once. By 
monitoring inclination and semi-major axis at the same time for 
a given planet, we can easily observe what influence the incli- 
nation has on the migration of the planet. The planets embedded 
in the isothermal disc are modelled with the e-potential, using 
r„„ = 0.8, while the planets embedded in the fully radiative disc 
are modelled using the cubic potential with r,,,, - 0.5. 
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Fig. 8. Density (in g/crn^) of a r - 0-slice through the isothermal 
H/r = 0.037 disc with embedded IQMEanh planets on fixed cir- 
cular inclined orbits with ; = 1.0° (top) and ; = 4.0° (bottom). 
The planet is at its lowest point in orbit (lower culmination) at 
the time of the snapshot. 

4.1. Isothermal disc 

Our disc only extends 7° above and below the discs midplane, so 
when simulating a planet with a higher inclination than that, it is 
not fully embedded in the disc any more. As the disc's material 
is concentrated mainly in the middle of the disc, the material 
high above or below the disc's midplane will only have a small 
effect on the planet anyway. Test calculations with an extended 6 
region of the disc (10° above and below the discs midplane) and 
an io - 8.0° 20 MEanh planet have shown that we do not need 
such a high extension above and below the disc. The aspect ratio 
for the disc is H/r - 0.037 in the isothermal case, as this aspect 
ratio corresponds to the aspect ratio of the fully radiative disc. 

In Fig. [TT] the evolution of inclination of a IQMEanh planet 
in an isothermal Hjr - 0.037 disc is displayed. After the plan- 
ets have reached their final mass after 10 orbits, the inclination 
begins to drop immediately. Up to / = 4.0° the planet's loss of 
inclination is increasing, while for even higher inclinations the 
damping of inclination is slowed down again (see Fig. [TJ. High 
incUned planets will first lose inclination at a quite slow rate, but 
as their inclination is damped the damping rate increases until 
the inclination reaches a; 4.0° and will slow down after that until 
the planets inclination is damped to zero. This results confirms 
our results for planets on fixed orbits. 

When the planet remains in the main body of the disc for 
i < Hjr (i in radians), the da mping of incUnation is exponen- 
tial: di/dt oc -/. By following iTana ka & WardI (|2004 the Hnear 
damping rate for a small planet mass and a small incUnation can 
be obtained: 




Fig. 9. Surface density for 20M Earth planets on fixed circular in- 
clined orbits embedded in a fully radiative disc with ; = 1.0° 
(top) and / = 4.0° (bottom). The planet is at its lowest point in 
orbit (lower culmination) at the time of the snapshot. 



with the characteristic time 



(12) 



In Tab.[T]we display t,„c for the isothermal and fully radiative 
simulations. We compare the linear estimate according to eqs.fTTI 
and[T2]with our full non-linear results from the fixed and moving 
planet simulations. 

For our isothermal Hjr - 0.037 simulations we estimate a 
linear damping time scale of r,„f = 14.2 orbits, which is about 
a third smaller than our obtained result T„,f = 23 for planets 
on moving orbits (see fits in Fig. fTTl i. The diff'erence between 
our actual fit and the linear value can be caused by the rela- 
tively small aspect ratio of the disc (///r = 0.037), because the 
linear theory is formally valid only for embedded planets with 
H » Rum- Below, we present additional result for anH/r - 0.05 





linear 3D 


fixed orbit 


moving orbit 


iso H/r = 0.037 


14.2 


27 


23 


iso H/r = 0.05 


47.25 




53 


fully radiative 


27.7 


20 


26.5 



di/dt 



1 



0.544 



(11) 



Table 1. The exponential damping rate T,„t. (in orbits), obtained 
from our isothermal and fully radiative simulations for three 
studied cases of planet s in discs c ompared with the linear rate, 
as given by[Tanaka &\\ ^(|2004 . Shown are results for planets 
on fixed orbits and for moving planets. 
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Fig. 10. Density (in g/cm^) of a r - ©-slice through the fully ra- 
diative disc with embedded IQMEanh planets on fixed circular 
incHned orbits with / = 1.0° (top) and ; = 4.0° (bottom). The 
planet is at its lowest point in orbit (lower culmination) at the 
time of the snapshot. 
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Fig. 11. Time evolution of inclination for lOMEanh planets with 
individual starting inclinations ranging from / = 1.0° to / - 8.0° 
in an isothermal H/r - 0.037 disc. The damping of inclination 
follows the trend of the calculated change of inclination for plan- 
ets on fixed orbits, see Fig.[T] The black dotted lines indicate our 
manual fitting with t,„c - 23 and in the io - 8.0° case an addi- 
tional fit provides for a di/dt cc fit. 



disc and find indeed better agreement of linear and fully non- 
linear results. However, our obtained r,„c = 23 for moving plan- 
ets is very close to the predicted damping time scale of t,„(. = 27 
for planets on fixed orbits, demonstrating the consistency of our 
results. 

The exponential damping la w should be only valid up to / ^ 
H/r, but lCresswell et all (l2007h stated that it is valid up to / 
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Fig. 12. Time evolution of semi-major axis for lOMeanh planets 
with individual starting inclinations ranging from ; = 1.0° to 
/ = 8.0° in an isothermal H/r - 0.037 disc. The dashed black 
line indicates a fit to the results with a slope of -3.45 10"^/orbit. 
The change of semi-major axis follows the trend of the change 
determined for planets on fixed orbits, see Fig. |3] 



2H/r, which in fact reflects our findings. The higher inclined 
planets do not lose inclination with an exponential rate at the 
beginning of the simulations, but with di/dt cc i^^. This is also 
indicated through the black dotted line in Fig. [TT]for the io = 
8.0° planet. As the planets lose their inclination, the damping 
becomes exponential again (at about / a; 3.0°). 

In Fig.[T2]the evolution of the semi-major axis in an isother- 
mal H/r - 0.037 disc is displayed. After a few orbits the semi- 
major axis decreases for all planets. During time the semi-major 
axis shrinks more and more, as you would expect from low mass 
planets in an isothermal disc. In the beginning the planets with 
io - 8.0° and io - 6.0° have a slower inward migration rate than 
the other planets, but as the inclination is damped in time, the in- 
ward migration settles for the same rate for all planets, although 
it takes about 80 orbits for the io = 8.0° planet. The initial higher 
inclination just delays the inward migration for a few planetary 
orbits, as could be expected from our results for planets on fixed 
orbits, see Fig.|3] 

The observed migration rate for all planets is constant af- 
ter the inclination has been damped. The black dashed line in 
Fig. [T2l indicates a linear fit for the evolution of the semi-major 
axis with a = -3.45 10"'* / orbit. This decrease of semi-major 
axis from planets on moving orbits is within 15% of the result 
for planets on fixed orbits in Fig. |3] The slight difference may 
be just due to the different setup, stationary versus moving and 
cubic versus epsilon potential. Additionally, the fit in Fig. [T2lis 
for a later evolutionary time, where the planet has already lost a 
few percent of its semi-major axis. 

In our past simulations we found a dependence of the migra- 
tion rate/torques due to the aspect ratio of the disc (IKlev et al.l 
120091: iBitsch & Klevl[20Toh . It is now logical to assume that the 
aspect ratio of a disc will change the damping of incUnation for 
inclined planets as well. In Fig.[T3]we display the evolution of 
the semi-major axis (top) and inclination (bottom) of planets em- 
bedded in an isothermal H/r - 0.05 disc. 

The inward migration of the planets in the H/r - 0.05 disc 
is slower compared to plan ets embedded in ?l H /r - 0.037 disc, 
which is in agreement with lTanaka et al.l (l2002l) . The inclination 
damping is slower in the H/r - 0.05 disc by a factor of 2 to 
3. Nevertheless, the final outcome is the same, an initially high 
inclination reduces the inward migration in the beginning, but as 
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Fig. 13. Time evolution of the semi-major axis (top) and inclina- 
tion (bottom) of a IQMEanh planet in an isothermal H/r - 0.05 
disc. The planets initial inclinations reach from /o = 1.0° to 
io = 8.0°. The dashed black line in the top panel indicates a 
linear fit to the results with a slope of -2 10""'/orbit. The black 
dotted lines in the evolution of inclination (bottom) indicate our 
manual fitting with t„,c = 53. 



inclination is damped, the inward migration settles for the same 
rate for all the inclinations (but with different rates for different 
aspect ratios). 

The theoretical, linear damping of inclination obtained by 
following Tanaka & Ward (2004) is r,,,,. = 47.25 orbits (eq.fT2li. 
which is slightly smaller than our fitted value of t,„c = 53 orbits 
(Fig- US bottom). The numerically obtained ratio of the migra- 
tion rates (53/23 ^ 2.3) differs from the theoretical linear ratio 
of (0.05/0.037)** a: 3.3. As the colder disc makes linear theory 
less applicable in these circumstances, we can infer that the re- 
sults obtained for the larger H/r - 0.05 disc are more accurate. 
The numerical exponential decay in the simulations with mov- 
ing planets matches our estimated rate for planets on fixed orbits 
quite well. The higher pressure inside the disc seems to have a 
smoothing effect on the damping, resulting in more accurate re- 
sults for discs with higher aspect ratio. 

For the moving planets, the measured migration rate (da/dt) 
for the H/r - 0.037 disc is approximately larger by a factor 
of (0.05/0.037)2 compared to the H/r = 0.05 disc (see dashed 
black line fit in Fig.[T3] top), which is to be expected in the linear 
case. 

In an isothermal disc, inclined planets lose inclination and 
semi-major axis immediately after they are released in the disc, 
as predicted by our calculations of the change of inclination and 
semi-major axis for planets on fixed orbits. The migration rate 
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Fig. 14. Time evolution of inclination for lOMEanh planets with 
individual starting inclinations ranging from ; = 1 .0° to ; = 4.0° 
in a fully radiative disc. The damping of inclination follows the 
trend of change of inclination for planets on fixed orbits, see 
Fig. [1] The black dotted lines indicate our manual fitting with 
Tine = 26.5 for planets with /q = 2.0°, 4.0°, 6.0° and /o = 8.0°, 
which seems very good in the beginning/middle phase of the in- 
clination damping. In the /q = 8.0° case an additional fit provides 
for a di/dt oc r'^ fit in the beginning of the evolution. 

is the same for all initial inclinations, a higher inclination just 
delays the inward migration for a few orbits. 

4.2. Fully radiative disc 

In the fully radiative regime we simulate planets to an inclination 
up to / = 8.0°. As always we use a 2D model in r-0 direction in 
the radiative equilibrium for the starting configuration of our 3D 
fully radiative dis c. This procedure is described in more detail in 
(iKlev et alj|2009t) . As in the isothermal case, the planet needs 10 
planetary orbits to reach its full mass. 

In Fig. [14] the change of inclination for moving planets in a 
fully radiative disc are displayed. The inclinations range from 
io - 1.0° to io = 8.0°. The planets inclination reduces as soon 
as the planet is released in the disc according to the rate found 
in Fig.[T] After about 160 orbits the inclination of all planets has 
reached zero, so that the planet is now orbiting in the equatorial 
plane. 

The estimated linear value for the damping of inclination in 
the radiative disc is T,„f = 27.76 orbits (eq.[T2li. which lies within 
about a few percent of our fitted value of t,„c = 26.5 orbits 
for moving planets. The difference with the calculated damping 
timescale r,„(. = 20 orbits for planets on fixed orbits is slightly 
larger in this case (see Tab. [TJ. But one should keep in mind as 
well that the exponential fit is valid for very low inclinations in 
the isothermal case as well, while it is only valid down to x 0.5° 
in the fully radiative scheme. 

On the other hand, the decay of inclination for initially high 
inclined planets fits even a little better than in the isothermal 
case, and the trend of an di/dt oc decay is clearly visible. 

In Fig. [15] the evolution of the semi-major axis of plan- 
ets with initial inclination is displayed. After the planets have 
reached their final mass at 10 planetary orbits, the initially lower 
inclined (/ up to 4.0°) planets start to migrate outwards as ex- 
pected from Fig. [3] The migration rates for planets on fixed or- 
bits also matches nicely with the determined migration rate of 
1.35 lO^-'/orbit for planets on moving orbits. The two higher in- 
clined planets first stay on their initial orbit before they start their 
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Fig. 15. Time evolution of semi-major axis for lOMEanh planets 
with individual starting inclinations ranging from / = 1.0° to 
i = 8.0° in a fully radiative disc. The change of semi-major axis 
follows the trend of change for planets on fixed orbits, see Fig.|3] 
The dashed black line in the top panel indicates a linear fit to the 
results with a slope of 1.35 10 Vorbit. 



outward migration. It seems that the planet has to lose a certain 
amount of inclination before it can start to migrate outwards. 
This result is agreement with our observations for planets on 
fixed inclined orbits. A higher inclination weakens the torque 
responsible for outward migration and if the inclination is too 
high, the outward migration stops. 

The planets seem to migrate outwards at nearly the same 
speed for all initial inclinations up to 4.0°. But before the plan- 
ets start to migrate outwards their initial inclination is damped by 
about 25% for all initial inclinations, which explains the similar 
migration speed for these planets, see Fig.|3] 

Low inclinations seem to have only little effect on the migra- 
tion of planets in the fully radiative scheme, but if the inclination 
of the planet is so high that its orbit is at some times high above 
the midplane of the disc, the disc becomes so thin that outward 
migration is not possible any more. 



5. Planets with different masses 

Small incUnations seem to have no big effect on migration for a 
2QMEarth planet on a circular orbit in the isothermal or fully ra- 
diative scheme. The planet's inclination is reduced in both ther- 
modynamic systems and it migrates inwards in the isothermal 
and outwards in the fully radiative scheme. In our previous work 
we found that plan ets up to x 33ME arih migrate outwards in a 
fully radiative disc (iKlev et alj|2009h . Planets with higher plan- 
etary masses migrate inwards, and considering our previous re- 
sults it is unlikely that inclination changes this. Nevertheless we 
here present simulations of planets with masses ranging from 20 
up to IQQMEartii in an isothermal H/r - 0.037 disc and in a fully 
radiative disc. We only present results for moving planets with 
our usual tool for increasing the planetary mass to its designated 
mass during the first 10 orbits. 

For all simulations we limit ourselves to the cases of /q = 
1.0° and io = 4.0° for all planetary masses, in the isother- 
mal H/r - 0.037 disc as well as in the fully radiative disc. 
Additionally, we display Jovian mass planets with inclinations 
up to 15° in the isothermal case. For the isothermal case, we 
again use the e-potential with a smoothing length of rj„, = 0.8, 
as it turned out to be more stable in the isothermal case than the 
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Fig. 16. Time evolution of the inclination for planets with io - 
1.0° (top) and io - 4.0° (bottom) for different planetary masses 
in an isothermal H/r = 0.037 disc. 



cubic potential, while we use the cubic r,,,, = 0.5 potential for 
the fully radiative simulations. 

In all simulations displayed in this section the planet is in- 
serted in the disc and starts to move immediately. This means 
that a massive planet will open a gap in time in the initially un- 
perturbed disc. If the planet were allowed to open a gap before it 
starts to move in the disc, the evolution of the planet, especially 
the evolution of inclination, would be slower since an open gap 
slows down migration and inclination damping. 

5.1. Isothermal disc 

In Fig.[T6]the evolution of inclination for planets with io = 1.0° 
(top) and /() - 4.0° (bottom) for different planetary masses in 
an isothermal H/r - 0.037 disc is displayed. The inclination 
drops immediately after the planet reaches its destined mass at 
nearly the same rate for all planetary masses in both cases, for 
the io = 1 .0° and io - 4.0° case. It seems that a higher planetary 
mass results in a slightly faster damping rate for the inclination. 
However, when the inclination reaches about 20% of the initial 
inclination, the 80 and lOOMEanh planets stop their fast inclina- 
tion damping for a few orbits. After about 150 orbits the inclina- 
tion of all planets is damped to zero, even for the initially more 
highly inclined planets. 

In Fig. [T7] we display the surface densities for planets with 
different masses for the io = 1.0° (top) and io = 4.0° (bottom) 
cases after 60 planetary orbits. The two surface density profiles 
are very similar, although the planets have a different starting in- 
clination io and a different inclination at the time of the snapshot. 
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Fig. 17. Surface density for planets on initial inclined, circular 
orbits with /q = 1.0° (top) and /q = 4.0° (bottom) after 60 plane- 
tary orbits. 



As the inclination of the /q — 4.0° and /q = 1 .0° planets evolves 
very similar (in the trend), it seems not surprising that the surface 
densities reflect this behaviour On the other hand, the surface 
density profiles after 60 orbits do not give a clear hint to explain 
the slowed down inclination damping for the high mass planets 
when they reach x 20% of the initial inclination. As expected, 
higher mass planets clear deeper and wider gaps inside the disc. 

The measured damping time scale r,„f for planets with dif- 
ferent planetary masses is displayed in Fig. [18] Starting from the 
smallest planet in our calculations, the lOMEwth, the damping 
time scale reduces until the planets mass is AQM Earth- For higher 
mass planets the damping time scale increases again. The in- 
crease for the 80 and lOOMgo,,/, planets should be handled with 
care, as the measurement for these planets is quite difficult due 
to the bump in the evolution of inclination. If we do not include 
these two planets in our thoughts, we can safely conclude, that 
the inclination damping increases with the planetary mass. The 
inclination of a more massive planet will be damped faster than 
the inclination of a smaller planet, even more so for radiative 
discs. 

In Fig.[T9]the evolution of the semi-major axis for the above 
mentioned planets is displayed. Interestingly it seems that the 
evolution of the semi-major axis does not depend on the incli- 
nation for circular orbits. All planets with equal masses follow 
the same evolution independent of the initial inclination. For the 
2QM Earth planct we found that the migration rate in the isother- 
mal case does not depend much on the inclination of the planet, 
see Fig. [3] For higher mass planets this feature seems to be true 
as well. Of course, the planets migrate inwards at a speed depen- 
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Fig. 18. Measured inclination damping time scale T,„f for plan- 
ets with different planetary masses in the isothermal and fully 
radiative regime. 
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Fig. 19. Time evolution of the semi-major axis for planets with 
(Q = 1 .0° (top) and /q = 4.0° (bottom) for diff'erent planetary 
masses in an isothermal HIr - 0.037 disc. 

dent on the planetary mass. A higher planetary mass results in a 
more rapid inward migration of the planet, as it is predicted by 
linear theory (Tanaka et al. 2002), although the theory formally 
applies only for an unperturbed disc. Our discs, however, show 
clear signs of perturbation in the surface density profile (Fig.fTTli. 

In Fig. |20] the evolution the inclination for planets with a 
Jovian mass (and for an additional planet with lOMEarth) in 
isothermal discs with different aspect ratios is displayed. The 
damping of inclination follows the predicted trend that the in- 
clination of planets embedded in discs with smaller aspect ratio 
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Fig. 20. Time evolution of inclination for Jovian mass planets in 
isothermal H/r = 0.037 and H/r - 0.05 discs. The evolution of 
a IQMEarih planct in an isothermal H/r - 0.037 disc is added 
for the inclination plot. The black dotted lines indicate a fit with 
Tine = 120 for a Jovian mass planet in an H/r - 0.037 disc and a 
fit with Tine = 145 for a Jovian mass planet in a H/r = 0.05 disc. 



(H/r - 0.037) is damped faster than the inclination of plan- 
ets embedded in discs with higher aspect ratio (H/r - 0.05), 
which was the result for small er mass planets in our pr evious 
isothermal simulations as well. Im arzan &Nelsonl(l2009h stated 
that inclination damping of higher mass planets (e.g. Jovian 
mass) is considerably faster than for low mass planets (e.g. 
2QM Earth)- Fig. [20l shows clearly the same result for high incli- 
nations (i'o = 15.0), as the 2QM Earth planet loses inclination at a 
much slower rate. 

However, for low incUnations (/ < 4.0°) the damping of in- 
clination for Jovian mass planets is considerably slower than for 
low mass planets. The damping time scales obtained by fitting 
our simulations are considerably higher for Jovian mass planets 
in isothermal discs, in the H/r - 0.037 disc Ti„e,ji,p = 120 > 
23 = Tine,2{)ME Well as in the H/r - 0.05 disc Ti„ejup - 145 > 
53 = Tine,2QME- The /() = 10.0° planet with IQMEarth reflects this 
behaviour very well. First the damping of inclination is slower 
than for a Jovian mass planet, then the damping of inclination 
increases and the planet ends up in the midplane of the disc. The 
damping of inclination strongly depends on the interactions be- 
tween disc and planet. As the Jovian mass planet opens a gap in 
the disc during the time of its evolution, the interactions between 
planet and disc are reduced, when the gap is opened. Therefore 
the damping of inclination of a Jovian mass planet is slowed 
down compared to a 20M£„,,/, planet, when the inclination is al- 
ready damped to about / ^ 4.0° after the planet has evolved for 
about 120 orbits. After the gap has opened, only little mass re- 
mains adjacent to the planet to provide continuing damping. In 
the beginning of the simulations the gap has not opened yet and 
inclination damping is faster. 

Overall, the damping of inclination starting from high in- 
clinations (io = 15°) is faster for planets with a Jovian mass 
compared to small mass planets {IQMEarth) as the damping 
time of inclination at high inclinations is much slower for low 
mass planets, which overcompensates the faster damping at 
low inclinations. This is in agreement with the statement of 
iMarzari & Nelsonl(l2009b . 

Inclination seems to have no visible effect on the migration 
of planets inside an isothermal disc, as long as the planet's orbit 
is circular. Eccentric orbits, on the other hand, tend to change the 
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Fig. 21. Time evolution of the inclination for planets with /q = 
1.0° (top) and /q = 4.0° (bottom) for different planetary masses 
in a fully radiative disc. 



migration rate of the planet. Later on, we also investigate planets 
on inclined and eccentric orbits at the same time. 



5.2. Fully radiative disc 

In Fig.|2T|the evolution of inclination for planets with /q = 1.0° 
(top) and /() - 4.0° (bottom) for different planetary masses in 
a fully radiative disc is displayed. The inclination is damped for 
all planetary masses, but a higher planetary mass favours a faster 
damping of inclination for both starting inclinations. After about 
150 orbits the inclination of all planets with both individual start- 
ing inclinations is damped to approximately zero. The damping 
rate of the inclination is about the same for the fully radiative 
and the isothermal H/r - 0.037 disc. 

In Fig.[T8]the measured damping time scale t,„c is displayed 
for isothermal and fully radiative simulations. In the fully ra- 
diative case, the damping time scale is reduced for increasing 
planetary masses. If the planet reaches a certain mass (about 
SOM^a,-,/,), the damping time scale is not reduced any more. As 
the planets start to open gaps in the disc, the damping of incli- 
nation is not further accelerated. It seems that gap opening plan- 
ets do not only migrate inwards at the same speed, independent 
of the planets mass, but also damp their inclination at the same 
speed. 

In Fig.|22]the evolution of the semi-major axis for planets 
with /'o = 1 .0° (top) and Iq - 4.0° (bottom) for different planetary 
masses in a fully radiativ e disc is displaye d. As expected from 
our previous simulations (iKlev et al.ll2009l) only planets with a 
planetary mass of Mp < 33MEarth migrate outwards, while heav- 
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Fig. 22. Time evolution of the semi-major axis for planets with 
io = 1.0° (top) and /o = 4.0° (bottom) for different planetary 
masses in a fully radiative disc. 

ier objects migrate inwards. The 35 and AQMEanh planets seem to 
migrate inwards just a little bit, before their migration is nearly 
stopped. The higher mass planets (50, 80 and lOOM^o,,/,) mi- 
grate inwards at nearly the same speed as the planets inside the 
isothermal Hjr - 0.037 disc. These observations are indepen- 
dent of the planets inclination. 

The outward migrating low mass planets travel outwards at 
a slightly faster speed in the /q = 1.0° case compared to the 
; = 4.0° simulations. This can be expected, as the calculated mi- 
gration rate for IOMecuiIi planets on fixed orbits is about a factor 
of 5 to 8 higher in the low inclined case. The higher migration 
rate for the lower inclined planets results in a faster outward mi- 
gration. Finding the observed difference in the semi-major axis 
being so small after about 150 orbits is a result of inclination 
damping. As the inclination is damped, the migration rate in- 
creases significantly until the inclination is ^ 2.0°, which takes 
only a few orbits starting from /q - 4.0°. In this little time frame, 
the planet is not able to migrate a large distance, so we only ob- 
serve a small difference in the evolution of the semi major axis. 

6. Inclined and eccentric planets 

Inclination did not change the direction of migration for plan- 
ets on circular orbits. Eccentricity on the other hand stopped the 
outward migration of low mass planets in fully radiative discs, 
as long as the eccentricity of the planet is higher than e - 0.02 
(iBitsch & Kiev 20 LQ). However, planets with lower eccentricity 
migrate outwards in the fully radiative scheme. If an inclined 
planet surrounds a star on an eccentric orbit, in what way will 



these two properties of the orbit influence each other? To find an 
answer, simulations of a low mass planet (lOMEanh) on an in- 
clined and eccentric orbit in an isothermal H/r - 0.037 and in a 
fully radiative disc are computed. 

6.1. Isothermal disc 

Our recent work with planets on eccentric orbits has shown that 
for isothermal simulations with moving planets on eccentric or- 
bits the cubic potential (for the planetary potential) should be 
replaced by the e-potential. In the following, isothermal simu- 
lations of planets on eccentric and inclined orbits we use the 
e-potential with a smoothing length of rj,„ = 0.8. 

In this work we limit ourselves to planets with an initial ec- 
centricity of eo = 0.1, eo = 0.2 and eo = 0.4 and with an initial 
inclination of /q = 2.0° and /q = 4.0°. In Fig.|23]the evolution 
of the semi-major axis (top), the eccentricity (middle) and incli- 
nation (bottom) for planets in an isothermal H/r - 0.037 disc is 
displayed. 

The planets migrate inwards (loss of semi-major axis) as ex- 
pected for low-mass planets in an isothermal disc. For planets 
with an initial inclination of /q - 4.0°, the planets migrate in- 
wards at a slower rate compared to the lower inclined cases. 
A lower initial eccentricity results in a faster inward migra- 
tion in the beginning of the simulation for all inclinations. This 
faster inward migration s eems to be in a direc t correlation with 
the planet's eccentricity dBitsch & Klevl[20Tol) . When the plan- 
ets eccentricity reaches about e ^ 0.1, the planet undergoes a 
short rapid inward migration and a more rapid loss of eccentric- 
ity, compared to the loss of eccentricity at different eccentric- 
ities. As soon as the eo = 0.2 planet's eccentricity is damped 
to about e ^ 0.1, the planet undergoes a rapid inward migra- 
tion leading to nearly the same evolution of semi-major axis as 
the eo = 0.1 planet from this time on, again for both inclina- 
tions. Interestingly, a lower starting inclination results in a faster 
inward migration, when the planets have an initial eccentric- 
ity, compared to planets on circular orbits, where the inclination 
seems to have no effect on the migration. A higher initial incli- 
nation (/ = 4.0°) seems to result in the opposite: a slower inward 
migration of the planet. We have observed a similar effect for 
eccentricities: planets with a lower initial eccentricity seem to 
migrate inwards at a slower rate compa red to their highly eccen- 
tric counterparts (IBitsch & Klevll20 10t). 

The initial eccentricity is damped for all these simulations 
with only two small diff'erences. Planets with an initially higher 
eccentricity need a longer time to be damped to zero and plan- 
ets with the same eccentricity but a higher inclination also need 
a longer time to damp their eccentricity. As the eccentricity is 
damped to lower values, at e a: 0.1 we observe a more rapid 
loss of eccentricity compared to the loss of eccentricity at other 
times. At this point the planet loses about 50% of its (actual) ec- 
centricity in the time of a few orbits. At this point in time we 
also observe a faster inward migration of the planet. This depen- 
dency of a rapid drop of eccentricity and a resulting decrease of 
semi-major axis in t he isothermal regime i s in agreement with 
the re sults found in ICresswell et al.l (l2007b and IBitsch & Klevl 
(120101) 

The evolution of the inclination for the low eccentric planet 
(eo = 0.1) follows nearly the evolution of a planet on circular 
orbit. The damping time is in fact a marginally longer for the ec- 
centric case, but besides this effect there is no difference. For the 
eo = 0.2 case the inclination starts to oscillate slightly in the be- 
ginning of the simulation until it is damped to about / a; 3.0° or 
/ a 1.5°. Then the inclination follows approximately the damp- 
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Fig. 23. Time evolution of semi-major axis (top), eccentricity 
(middle) and inclination (bottom) for planets {2QM Earth) with 
different eccentricity and inclination in an isothermal H/r - 
0.037 disc. 
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Fig. 24. Time evolution of semi-major axis (top), eccentricity 
(middle) and inclination (bottom) for planets (lOMEanh) with 
different eccentricity and inclination in a fully radiative disc. 



ing rate of a planet on a circular orbit (see Fig. [TTl until it is 
damped to about zero. The high eccentric case (eo = 0.4) shows 
even more and stronger oscillations in the inclination. In the be- 
ginning the planets inclination is pumped up to about ^ 5.0°, 
which is 25% more than the starting inclination. In time, the am- 
plitude of the oscillations becomes smaller and the frequency 
higher while the overall in clination dissipa t es in ti me. These os- 
cillations are observed in ICresswell et al.l (l2007l) as well. The 
time until the inclination reaches zero is about 4 times as large 
as in the zero eccentricity case. 

By comparing the migration rate of the inclined and eccen- 
tric planet with only an eccentric planet (migrating in midplane 
with a zero inclination), we observe nearly the same migration 



rate for both planets. Due to these results, it seems that incli- 
nation does not influence the migration in an isothermal disc at 
aU. 



6.2. Fully radiative disc 

For the fully radiative simulations we rely on our cubic rj,,, - 0.5 
potential again. For comparison, the starting eccentricity and in- 
clination for the fully radiative simulations match those of the 
isothermal simulations. In Fig. |24] the evolution of the semi- 
major axis (top), the eccentricity (middle) and inclination (bot- 
tom) for planets in a fully radiative disc are displayed. 

In the beginning of the simulations all planets migrate in- 
wards in the same way as in the isothermal simulations, but in 
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contrast to the isothermal simulations the planets reverse their 
inward migration and migrate outwards after about x 100 or- 
bits in the low eccentric case (eo < 0.2). The high eccentric case 
(eo = 0.4) takes a much longer time to reverse its inward migra- 
tion. A higher initial eccentricity results in a later outward migra- 
tion of the planet. The evolution of the semi-major axis contains 
at some point a quite fast loss of semi-major axis for the initial 
low eccentric planets. This loss of semi-major axis correlates to 
the eccentricity. As soon as the eccentricity is damped to x 0.1, 
the fast loss of the semi-major axis sets in. The inclination of the 
planet seems to play no significant role in the migration process 
of these simulations, as for the eo = 0.1 case the higher inclined 
planet (/o = 4.0°) migrates outward first, while for the eo = 0.2 
case the lower inclined planet (/q = 2.0°) migrates outward first. 

The initial eccentricity is damped for all planets, but for plan- 
ets with the same starting eccentricity, the time needed to damp 
the eccentricity is longer for the higher inclined planets, see the 
middle of Fig. |24] As the eccentricity is damped to » 0.1, it is 
damped with an exponential decay, which seems to result in a 
faster loss of semi-major axis at the same time in the evolution 
of the planet. The damping for the initial low eccentric planets 
(e < 0.2) takes about 120 planetary orbits, while a higher initial 
eccentricity (eo = 0.4) takes about three times longer 

The inclination damping follows in principle the evolution of 
inclination of the isothermal simulations. As soon as the planet 
is released in the disc, the inclination is damped as long as the 
eccentricity is very low (eo = 0.1). Higher eccentricities seem to 
affect the inclination in such a way that the inclination begins to 
oscillate. This oscillation is already visible in the eo = 0.2 case, 
but becomes much stronger for even higher eccentricities. The 
oscillations begin to dissipate as soon as the planets eccentricity 
is damped to about e < 0.1. During time these oscillations lose 
amplitude and gain frequency (as we can see from the eo = 0.4 
simulation). A higher eccentricity seems to cause a higher and 
more frequent oscillat ion of the inclination. These oscillations 
have been observed in ICresswell et al] (l2007l) as well. 

7. Summary and Conclusions 

We performed full 3D radiation hydrodynamical simulations of 
accretion discs with embedded planets of different masses on in- 
clined orbits. In a first sequence of simulations we have analysed 
in detail the dynamics of a planet with a given mass of lOMEanh 
for the isothermal as well as the fully radiative case. For plan- 
ets on fixed orbits we calculate the expected change of inclina- 
tion and semi-major axis of an inclined planet on a circular orbit 
embedded in the disc, which we confirm subsequently through 
simulations of moving planets. For the isothermal {H/r = 0.037 
and Hjr - 0.05) as well as for the fully radiative simulations 
the inclination is damped in time. The damping time scale in 
isothermal discs depends on the disc's thickness {Hjr), as the in- 
clination is damped faster for discs with lower aspect ratio than 
in discs with higher aspect ratio. A smaller Hfr also results in a 
faster inward migration compared to a higher H/r . Om s imula- 
tion results agree in this respect with lTanaka et al.l(l2002l) . 

For planets on inclined orbits we confirm the outward migra- 
tion. Our results for planets on fixed orbits agree very well with 
the simulations of planets moving inside the disc. We find that 
outward migration interestingly occurs over a relatively large 
range of inclination, up to about (/ = 4.5°), with a faster outward 
migration for lower inclined planets. For higher inclined planets 
the migration seems to be stalled in general. In our simulations, 
the damping of inclination is slower in the fully radiative disc 
(Tine = 26.5) compared to the isothermal disc having the same 



scale height (T,„e = 23), a difference that may be attributed to the 
different sound speeds, isothermal versus adiabatic. 

As recent studies showed, planets in a fully radiative disc 
only migrate outwards if they do not exceed a certain critical 
mass (iKlev et alJ l2009h . beyond which gap formation sets it. 
This outward migration does crucially depended on the shape of 
the orb i t of the planet in t he disc, as discussed in lBitsch & KlevI 
(l2010l) : iKlev et all (l2009h . Here, we extended those studies to 
inclined planets and found that planets on circular and inclined 
orbits migrate inwards in the isothermal case, and migrate out- 
wards (if the planetary mass is low enough and the inclination 
of the planet is below the threshold of / w 4.5°) in the fully ra- 
diative case. For higher mass planets, our previous results have 
been confirmed in the expected way. Planets with a mass up to 
~ 33MEarth migrate outwards and planets with higher masses 
migrate inwards in the fully radiative scheme. In the isother- 
mal case, planets m igrate inwards as predicted by linear theory 
(iTanaka et al.ll2002l) . When evolving planets with higher plan- 
etary masses, the simulations started from unperturbed discs, 
which changes the damping rate of inclination compared to plan- 
ets starting in perturbed discs (by the presence of the embedded 
planet). This was a lso stated for inclined hig h mass planets in 
isothermal discs bylM arzan &Nelsonl(l2009r After the inclina- 
tion is damped in time, the direction of migration only depends 
on the thermodynamics of the disc and on the planetary mass. 

For eccentric and inclined low mass planets we observe a 
quite different behaviour Eccentricity seems to slow down the 
inward migration of planets in the isothermal regime and pre- 
vents outward migration in the fully radiative scheme (for low 
mass planets which would normally migrate outward). Since ec- 
centricity is damped with time, this leads to slightly faster inward 
migration in the isothermal regime and to a reverse of the migra- 
tion direction in the fully radiative scheme. The origin of this 
outward migration is created by the delicate density structure 
near the planet, which is destroyed b y even a very low amount 
of eccentricity (iBitsch & Klevll20Tol) . A high initial eccentricity 
gives rise to oscillations in the inclination, which become weaker 
and more frequent in time, as the eccentricity is damped. These 
oscillations seem to be responsible for a slower damping of in- 
clination compared to low eccentric planets. 

The results by iMarzari & Nelsoni (l2009b for Jovian mass 
planets indicated that inclination and eccentricity is damped in a 
very short time compared to low mass planets. We can confirm 
their results for inclination damping in discs with two different 
aspect ratios, however this is only true for the damping of high 
inclinations. For low inclinations (; < 10°) the damping of incli- 
nation for low mass planets is faster than for high mass planets. 
However, when considering the total time needed to damp the 
inclination of a planet starting from /o - 15.0°, the damping is 
much faster for high mas s planets, which confirms the results by 
Im arzan & Nelsoni (l2009h . A Jovian mass planet with high incli- 
nation evolves more rapidly as its interactions with the disc is 
stronger if no gap has opened yet. 

Our results indicate that the inclination as well as the eccen- 
tricity of single planets are damped by the disc. At the end of our 
simulations the planets end up in midplane, meaning they have 
no remaining inclination. While this is in agreement with the 
very flat solar system, it seems to be in contradiction to the ob- 
served highly inclined exoplanets (iTriaud et al.ll2010l) . However, 
it is known that higher eccentricities as well as inclinations can 
be produced via planet-planet interaction. For a pair of planets 
embedded in discs convergent migration processes can excite ec- 
centricity toJiigh_valuesdepending on the damping action of 
the disc (ICrida et al.ll2008h . Stronger scattering events can be 
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induced in case of a multipl e planet system still embedded in 
the disc (iMarzari et al.ll2010l) . During the subsequent longterm 
evolution these initially excited system may then evolve into 
confi gurations with very h i gh eccentricity as well as inclina- 
tion jChatteriee et al.l 120081: iMatsumura et al.ll201d) . Thus, the 
combination of initial planet-disc interaction with subsequent 
scattering processes and tidal interaction with the central star 
is one pathway t owards the observed strongly tilted systems 
( Winn et al1l2010l). O n the other hand, systems such as Kepler-9 



Holman et alJl2010h seem to show clear evidence for continued 



migration towards the star. 
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